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Boson clouds around black holes exhibit interesting physical phenomena through the Penrose 
process of superradiance, leading to black hole spin-down. Axionic clouds are of particular interest, 
since the axion Compton wavelength could be comparable to the Schwarzschild radius, leading to the 
formation of "gravitational atoms" with a black hole nucleus. These clouds collapse under certain 
conditions, leading to a "Bosenova" . We model the dynamics of such unstable boson clouds by a 
simple cellular automaton and show that it exhibits self-organized criticality. Our results suggest 
that the evolution through the black hole Regge plane is due to self-organized criticality. 
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I. INTRODUCTION 

Modern particle detectors are typically associated with 
very high energies, like the current design of the LHC 
with 8 TeV center of mass energy. The need for having 
high energies is obvious if the particles to be discovered 
— like the Higgs — have a substantial rest mass or, con- 
versely, a tiny Compton wavelength. 

Interestingly, for (QCD-)axions [lJ-Q the situation may 
be reversed: the rest mass of axions can be tiny, while 
their Compton wavelength can be substantial. More pre- 
cisely, the axion Compton wavelength can be of the same 
order of magnitude as the Schwarzschild radius of a stel- 
lar mass black hole. Thus, at least in principle, black 
holes can serve as particle detectors for axionsJ4, 5]. This 
idea has engendered a lot of recent interest [6-17]. 

The present work is based on the working assump- 
tion that suitable axions exist in our Universe (we shall 
be more precise what "suitable" means in the body of 
the paper). In that case an axion cloud forms around 
a rotating black hole, and the combined system black 
hole/ axion cloud can be thought of as a gigantic "grav- 
itational atom", with the black hole playing the role of 
the nucleus and the axion cloud playing the role of the 
electrons. 

As explained in [5j , the dynamics of this "gravitational 
atom" is governed by the Penrose process of superradi- 
ance [lg|. This process leads to exponential growth of 
the occupation numbers of the "atomic" levels, i.e., the 
formation of an axion Bose-Einstein condensate. The 
energy required for the formation of this condensate is 
extracted from the black hole, which thereby spins down 
and also loses some of its mass M . In fact, the black hole 
parameters (M, a) move on a Regge trajectory, where 
a £ [0, 1] is the Kerr parameter. Eventually, when the at- 
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tractive axion self-interactions become stronger than the 
gravitational binding energy, the axion cloud collapses — 
a "Bosenova" occurs. 

The aim of our present work is to map relevant as- 
pects of the dynamics of the black hole/axion cloud sys- 
tem to a simple cellular automaton, and to show that it 
exhibits self-organized criticality. If this conclusion holds 
also for more realistic cellular automata it would mean 
that self-organized criticality governs an essential part of 
the dynamics of gravitational atoms. 

This paper is organized as follows. In section |TT] we 
review some preliminaries, namely crucial properties of 
self-organized criticality and the black hole/axion cloud 
system. In section IIIII we present the working assump- 
tions of relevance for our cellular automaton and argue to 
what extent they are plausible. In section HVl we present 
results based upon computer simulations implementing 
this specific cellular automaton. In section [V] we summa- 
rize and point to open issues. 



II. PRELIMINARIES 

A. Self-organized criticality 

The concept of Self Organized Criticality (SOC) was 
proposed as an organizing principle of Nature, as it 
was realized that numerous spatial extended systems, in 
physics, biology, social media and economics, exhibit a 
number of properties which may be shortly character- 
ized as flicker noise for the temporal evolution and self- 
similar (fractal) behavior for the spatial evolution [19j . 
Ever since this seminal paper, the analytical and numer- 
ical methods developed to be able to tackle SOC in var- 
ious contexts have led to progress in understanding the 
nonlinear dynamics of complex, interacting systems. It is 
expected that SOC develops in slowly driven interaction 
dominated threshold systems and it has been shown that 
this state is an attractor of the dynamics of such systems 
(see for instance [2(| and references therein). 



Due to the complex interactions between its con- 
stituents, after a transitory phase in its evolution, the 
system exhibits what may be called a phase transition to 
a state characterized by the absence of a critical length 
scale, i.e., in this state all lengthscales are relevant. 

The archetypal SOC system is driven by an external 
force with a very large characteristic timescale and the 
system evolution is determined by local rules, i.e., only 
near- neighbor interactions are considered [20, [2l| . One 
important point to make is that this external driver is not 
artificial/conscious but is intrinsic to the physics in the 
system. This property of SOC is unlike the phase transi- 
tions in an Ising type system, where the phase transition 
is acquired by a conscious tuning of parameters to a crit- 
ical value. 

Because this line of research is still in its infancy and 
due to the associated analytical complexity, analysis of 
such systems relies heavily on numerical methods and 
so they have become intertwined in an inseparable way. 
It is difficult to enumerate properties of SOC without 
relating them to the way they are simulated in computer 
experiments. 

Systems which exhibit SOC are discrete (in general 
modelled by lattice models). The interaction between 
constituents are usually nearest neighbor-like and gov- 
erned by local rules and are modelled as interaction be- 
tween sites in a lattice, such that if one site receives an 
information, only a few other selected sites have access 
to this information. The entire system is subjected to an 
infinitely slow external drive and due to this drive and 
the complex interaction, the onset of SOC usually oc- 
curs when all the sites in the lattice host a near- threshold 
value of some parameter. If the threshold is reached, the 
system becomes subjected to different rules of evolution 
(critical flow) which allow information to be transmitted 
not only locally, but across all lengthscales. If one site has 
reached threshold and sends information to the next site, 
this next site might also reach threshold and so on. The 
avalanches of information from one seed site to a final site 
(i.e. where the threshold condition is no longer satisfied) 
form an event. The size of an event is the total number 
of avalanches occurring from one seed site. Natural sys- 
tems that reach SOC maintain it for an indefinite time 
as long as the external driver keeps acting. This equi- 
librium state is maintained by a balance between driving 
and dissipation, which occurs at the boundaries or in the 
bulk. 

The way that SOC was historically "detected" in Na- 
ture was due to the properties exhibited by observables of 
the system. These observables show a powerlaw distribu- 
tion of values. It was shown [2fJ,[21| that, for a system in 
SOC, the distribution function D(s) of some parameter 
s is given by 



D(s) 



y \bL") 



(1) 



where Q{x) is the universal scaling function for a given 
universality class of systems, a is a system specific con- 



stant, r is the critical exponent, b is a system dependent 
amplitude, L is the system spatial size and a is deter- 
mined by the physical dimension of the parameter s. A 
very useful feature of SOC is that if the distribution of 
another observable, say y, is considered, then D(y) will 
be equivalent to Eq. ([1}, with s replaced by y and dif- 
ferent parameters a, r, b and a, but the same scaling 
function Q. 

This is conceptually a powerful tool for observational 
purposes, in cases where it proves easy to record -D(s), 
but difficult to record a perhaps more interesting distri- 
bution D(y). 

An example of a system in SOC and a discussion of the 
properties emphasized above is given in the Appendix. 



B. Massive boson instabilities near a rotating black 

hole 

It has been shown 0, LL2|, [22| that the presence of 
bosons near a rotating black hole (BH) will lead to the 
formation of a boson cloud (BC) around the BH if a 
series of conditions are met. This BC-BH system is anal- 
ogous to an atom. The cloud is being continuously fed 
with bosons due to the superradiance phenomenon. At 
some point, when the number of bosons in the cloud is 
too large, the nonlinearities in the BC lead to an unsta- 
ble behavior and the cloud goes through a Bosenova-type 
event. 

The conditions for this chain of events to occur re- 
fer to a set of parameters characterizing both the boson 
and the BH, {M, w/j, (J-^r, rn}, where M and u>h are the 
mass of the BH and the angular velocity at the hori- 
zon, respectively, and /i, lor and m are the mass of the 
boson, angular velocity of the associated wavepacket and 
"magnetic" quantum number of the boson. The first con- 
dition requires that a boson exists such that its associ- 
ated Compton wavelength, Ah, is approximately equal to 
the characteristic size of the BH, R g (the Schwarzschild 
radius). It was shown that the QCD axion can satisfy 
this condition [in e.g. [5]]. Superradiance, i.e., contin- 
uous feeding of the cloud with "free" bosons, the phe- 
nomenon analogous to the Penrose process for fields, 
occurs if the superradiance condition is met, namely if 
cj e [0,u>hTn]. Once superradiance sets in, the topology 
of the growing cloud is given by a frequency ujr, which 
satisfies the superradiance condition and is the frequency 
of hydrogenic-like bound levels. Growth of the BC, in- 
stability and subsequent Bosenova occur if the frequency 
associated to the boson as it interacts with the BH has 
a very small imaginary part, T g , such that the frequency 
describing the wavepacket in the potential well of the BH 
is lu = ujr + iTg. The parameter r g , giving the growth of 
the BC, is a function of the set of parameters mentioned 
above and it exhibits a maximum for specific values of 
these parameters Q. 

We are now in the position to define the "suitable 
axion" mentioned in the Introduction. This axion is a 



boson which possesses physical characteristics such that 
Af, ~ R g , lur £ [O^uj^m] and, in the potential well of 
a BH, exhibits a growth rate T g which is very small, 
ljb. ^S> T g . The growth rate also depends on details of 
the bound system, such as the hydrogenic-like quantum 
numbers, taken here such that T g has its maximum pos- 
sible value. We assume that this boson exists and give 
a general outline of the dynamics of the BC-BH system. 
Superradiance sets in and the number TV of such bosons 
in the vicinity of the BH is amplified, in the linear limit, 
as 



dN 



T g N. 



(2) 



However, during this amplification more complex pro- 
cesses occur. Most importantly, the nonlinearities start 
increasing and the bosons begin to interact with each 
other. The end result is that the number of axions 
present in the cloud is reduced, with a rate T^ (d from 
dissipation), where Td <^T g . For our purposes, Eq. ([2]) 
can be naively modified as 



dN 
~dt 



TgN 



T,N. 



(3) 



Even in this very simple description it is possible to 
recover some of the qualitative results derived rigorously 
in [5j. Since dissipation does not equilibrate growth, the 
number N of bosons grows until the back-reaction on 
the frequency ujr is such that the shape of the cloud 
departs from hydrogenic (nonlinear limit). This has been 
shown [3, LL2| to occur when 
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(4) 



where Mbc is the mass of the boson cloud. 

The cloud grows further, but in a nonlinear regime, un- 
til the nonlinearities in the boson self interaction poten- 
tial are too large and the cloud is no longer stable. A large 
percentage of the cloud is supposed to collapse [9, LL2[ 
(critical limit). This collapse occurs on a timescale of r g 
of the BH, meaning that from the moment the critical 
limit has been reached, it takes r g units of time until a 
significant part of the cloud collapses into the BH. 

The subtlety here is that this process occurs following 
an avalanche of information spreading on all lengthscales 
accessible to the system (i.e. with an upper bound given 
by the system size). 

If the superradiance condition is still met, the cloud re- 
grows within the same parameter space {M, u>h, U, to, m}. 
In a plot of the BH spin vs. M\x (a Regge plot, @, their 
Fig. 3), the subsequent collapses and regrowths, with the 
same values of the parameters that initially satisfied the 
superradiance condition, resemble an almost straight line 
on which the systems wanders (Regge trajectories). It 
stays on this line for tens to hundreds of e-fold times of 



boson cloud growth, i.e. for tens to hundreds of Boseno- 
vas. When the superradiance condition is no longer met 
with the initial parameters, the system makes a transi- 
tion to another set of parameters (i.e., a new line on the 
Regge plot) which are again superradiant and analogous 
dynamics occurs. 



III. CELLULAR AUTOMATON FOR BOSON 
CLOUD AROUND BLACK HOLES 

We want to map the dynamics of BC-BH system into 
SOC. As was mentioned previously, the concept of SOC 
and the simulation of SOC are intimately connected and 
the only contained way to map the BC-BH dynamics to 
SOC is to simultaneously discuss the computer simula- 
tion associated with this analysis. SOC simulations are 
generally done using Monte Carlo methods on discrete 
grids, and we will focus on a Cellular Automata (CA) 
approach. The CA is a "machinery" which stores values 
of a parameter on a discrete grid and evolves these values 
in time and space according to simplistic rules designed 
to mimic a real phenomenon [23] . 

The CA used in the remainder of this paper is 2- 
dimensional cartesian grid, of size L x x L y . While this 
undoubtedly simplifies the numerical effort, there is also 
strong physical argumentation in favor of considering a 
2D grid (instead of a 3D grid). The BH has a signifi- 
cant spin, the quantum numbers of the atom are fixed 
(by assumption of the "suitable axion") and the orbits 
are Keplerian, which means that angular momentum is 
conserved thus leading to planar motion. 

Three questions then arise: 1) what is the quantity 
which will be placed on the grid, 2) whether it is physi- 
cally acceptable to discretize that quantity in this manner 
and 3) what are the evolution rules? We proceed by an- 
swering the first two questions here, while the third one 
will be the subject of Section HVl 

Our simulation procedure is based on the following rea- 
soning: if the number of axions in the cloud can be com- 
puted by some means [as in Q, their Eq. (28)], a local 
quantity can be defined by dividing the number of axions 
by the area of the cloud. Because the area of the cloud is 
considered to be constant, simulating the evolution of the 
density means simulating the evolution of the number of 
axions, and this is why we will refer to the quantity in 
the grid as "the number of bosons" . 

This basically means that the N bosons in the L x L 
cloud will be accommodated into a L x L grid, where we 
will from now on denote L x = L u = L. So it follows that 
one cell in the grid has A' = N/L 2 bosons and N will be 
the parameter stored in the CA. 

To answer the more important question if this dis- 
cretization is physically acceptable, consider that the 
number of bosons, N, or more precisely the N very 
poorly localized wavefunctions allowed by the superra- 
diance, may be replaced by wavefunctions which peak 
sharply, i.e. are more localized for simulation purposes 



(Fig. [TJ. This "effective" system has the same informa- 
tion content from the point of view of an outside viewer 
looking at gravitational radiation following collapse. One 
will see in the following algorithm/simulations that the 
information released following collapse does not depend 
on the actual position of where the information was re- 
leased. Thus, from the point of view of an external ob- 
server equipped only to record information on gravita- 
tional radiation, it is not important if an axion left the 
cloud by withdrawing its wavefunction from the entire 
cloud or only from a more localized spatial extent as in 
our simulation. 



TABLE I: Map of BC-BH dynamics to SOC. 




FIG. 1: Schematic representation of the distribution of infor- 
mation about the density in the cloud. Black lines: N poorly 
localized bosons. Blue lines: N "effective" bosons, better lo- 
calized for simulation purpose. 

We address now how the temporal evolution of the 
BC-BH system looks like in this framework. The linear 
dynamics of the BC-BH are implemented and are the 
evolution rules of the CA. N is increased at randomly 
chosen positions at each timestep in order to simulate 
superradiance and Nij(k) is the value of the variable at 
the grid-cell labelled by spatial indices i and j, at time 
step k, as explained in detail in Section llVl We record the 
value of a given parameter s characterizing the system 
and check if the distribution function of this parameter, 
D(s), behaves as for a SOC system. 

We summarize our assumptions 

• There is only one populated energy level. This is 
equivalent to saying that Eq. ([3]) expresses, in a 
linear approximation, the dynamics of the cloud 
itself. Such an assumption is justified, because the 
growth rate of the next available energy level is very 
small compared to T g . 

• The dissipation rate T^ is constant as the cloud den- 
sity grows. By assumption, dissipation is negligible 
compared to growth, so even if Td does change, it 
will not change as much as to affect the dynam- 
ics @. 

• The nonlinear and critical limit are the same and 
occur when the parameter e reaches its critical 



SOC/CA framework 


BC-BH system 


discrete space 


2D "cloud" 


local interaction 


one cloud cell interacts with 3 
other cloud cells only at criti- 
cality 


infinitely slow external drive 


bosons are added to the cloud 
with a very slow rate, T g 


SOC occurs as a result of 
threshold dynamics 


when the critical parameter 
e = 10 -4 is reached, nonlinear- 
ities become important 


dissipation 


given by the very slow rate Td 


observables 


distribution of the number of 
avalanches needed to relax one 
perturbation 



value 10 4 , Eq. (jU). This is supported by simu- 
lations of collapse [12[ • 

• the size of the cloud, L, is constant and is calculated 
based on the linear approximation of hydrogenic 
states. 

The mapping between BC-BH and SOC/CA proposed 
in this paper is summarized in Table HI 



IV. RESULTS 

We address now the third question posed in the pre- 
vious Section by formulating the evolution rules for our 
CA. Let L be the cloud radius, i.e., the r c quantity from 
Eq. (11) of @ 



(n + l + 1) 1 
{nR g f 



R„ 



(5) 



with n = 0,1 = 1, [iR g = 0.42 and fi = 3 • 10" n eF. n 
and I are quantum numbers associated with the hydro- 
genic state of the cloud and these particular values assure 
maximum growth. In this case, the constant [24J most 
effective instability rate occurs when [5>] 



r 9 =1.5-10- 7 rg- 1 . (6) 

The solution of the differential equation ([2]) is given by 



N(t) = N t e r ^, 



(7) 



where 2Vj is the initial value of the variable N. For sim- 
ulation purposes, temporal evolution is recorded in dis- 
crete timesteps, At. The solution at timestep k may be 
rewritten as 



N(k ■ At) = N % e 



r a kAt 



(8) 



If At is taken to be of the order of one gravitational 
radius and the notation T g = T g At is used, the equation 
for the evolution of the system in natural time units is 
formally written as 



r„fc 



N(k) = AV 9 



(9) 



The number of bosons added to the cloud due to su- 
perradiance in one natural timestep is given by 

AA = N(k + 1) - N(k) = N e fsk (e r * - l) . (10) 

Since for the QCD axion e kTs = 1 and e r » — 1 = 
1.5 • 10~ 7 , the dynamics equation used as a basis for the 
simulation is 



N iij (k + l) = N id {k)-N , 



(16) 



for m randomly chosen sites, with m <C n. For n = 100 
we consider m = 1. For future reference we note that in 
numerical simulations the important quantity is actually 
the ratio N c /N 



^ = 6.67 X aT^O 4 . 

An 



(17) 



To summarize, for a QCD axion with fiR g = 0.42, 
initial cloud population with x — 95% and after rescaling 
the values of A c and An. so that No is of order unity (for 
practical reasons) we get the following parameters used 
in the simulation 



N(k + 1) = N(k) + 1.5Ni ■ 10" 



(11) 



It is useful to restate this equality in terms of the critical 
number of bosons, N c . We will do this by assuming that 
the initial number of bosons is some fraction x £ (0, 1) 
of the critical number, A,- = xN r 



N(k + 1) = N(k) + 1.5xN c ■ 10" 



(12) 



id 



A c is easily computable as A c = Mbh^I V = 2.4 x 10 
(all quantities are calculated for mass expressed in eV). 
So the equation for the temporal evolution of the total 
number of bosons in the cloud can be written as 



A(fc + 1) = A(fc) + AA, 



(13) 



where AN — 3.6a; • 10 9 . This is to be simulated on a 
100 x 100 grid with the help of the parameter N. One 
may immediately write 



N iJ (k + l)=N lJ (k) + AN, 



(14) 



for all sites labelled i,j in the grid and at each timestep 
k. This is not general enough because so far we have no 
reason to believe that the spatial growth of the number 
of bosons is homogeneous. To circumvent this problem, 
the same total amount AN will be added to the cloud, 
but to a limited number of randomly chosen sites, i.e. 
the update rule for the CA will read 



N l , J (k + l)=N h3 {k)+N , 



(15) 



for n randomly chosen sites, with An. = AN/n. If n = 
100 then An. = 3.6a; • 10 7 . For the discretized case the 
critical parameter is N c — 2.4 • 10 12 . Dissipation must 
also be taken into account, but at a much smaller rate. 
This is modelled as 



iV c = 280000, N =4: 
n = 100, m = 1. 



(18) 



In words, at each time step k, a number 100 units are 
added to the cloud and iV grows locally. The local en- 
hancement of axion number is done randomly, i.e., 100 
randomly chosen grids will receive one unit. At the same 
timestep, 1 unit dissipates from a randomly chosen site 
and it is assumed that this unit escapes to infinity (it is 
not fed to the BH) . Adding new axions at randomly cho- 
sen positions is a nontrivial assumption, but it seems well 
justified based on two considerations: the observations of 
laboratory Bosenovas which emphasize the stochastic na- 
ture of this process 24| and the previous statement that, 
until the critical limit is reached, the system evolves in a 
linear regime, Eq. ([3]). 

At each timestep k all the grids are verified to see 
whether the critical density N c is reached. Once the crit- 
icality condition is met, the gravitational potential is no 
longer important and transfer of information is due to the 
nonlinear interaction of the bosons. The next question 
is "how many neighbors does the site interact with?" so 
that information of collapse could be transported to that 
number of neighbors. We will assume it interacts with 
3 sites, namely those three which are nearest neighbors 
and closer to the BH than the current site. The critical 
condition and the critical flow are implemented as 



ii Nj(k) > A c ,then 

Ni+i,j(k) -► N i+1 j(k) + N a ; 
iV"i+ij_i(fc) -> N i+ ij^i(k) + No; 



(19) 



N, 



i+lj'+l 



(fc) -> N 



i+lj'+l 



(A) + ^o; 



The result of this simulation is the distribution D(s) 
of event sizes (Fig. ^. This is exactly how an archetyp- 
ical SOC system behaves [19| and it is clear that after a 



growth period Self Organized Criticality is reached and 
the system stays in this state. This is a first important 
result. 
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FIG. 2: The distribution of event size, D(s) for a QCD-axion 
cloud in superradiance around a BH. The flow at criticality is 
given by the rules (I19|) . The cloud was initialized at 95% of 
its critical boson number N c . 



To test whether or not the fact that the system reaches 
SOC is rule-dependent, the system evolution is studied 
for different rules at criticality (different critical flows). 
The difference is that apart from the three units allowed 
to move to other sites, an additional number of boson 
units completely leave the cloud. This was done for 3 
units, 6 units [Eq. fl20|), Fig. |3] and 12 units. For the 
case of six units, this is written and implemented as 



if N i; j(k) > N c ,then 

N id {k)^N id {k)-9N ] 

N i+1J (k) ^ N i+ltj (k) + N ] 



(20) 



Ni+u-^k) -+-/V, 



Ni+i, J+ i{k) -+ A^+ij+i 



.. . i j_i(fc) + N ; 

(k) + N ; 




50000 100000 150000 200000 250000 300000 350000 400000 450000 50000( 

FIG. 3: The distribution of event size, D(s) for a QCD-axion 
cloud in superradiance around a BH. The flow at criticality 
is given by rules (|20l) . The cloud was initialized at 95% of its 
critical boson number. 



for Mbh = 0.9M. All the parameters characteristic to 
the system, and thus to the simulation, change accord- 
ingly. Simulations including the change in superradiance 
parameters were performed for all 4 sets of rules at crit- 
icality. The same qualitative results were obtained in all 
simulations: there is SOC on both sides of the transi- 
tion, the transition is sharp and the ratio of the number 
of events prior to transition and following transition is 
approximately 2. 

When exiting the superradiance condition for the pa- 
rameter set governed by M, superradiance for the param- 
eter set governed by — 0.9M instantly begins according 
to our rules. The results of this is that N, which has 
a value around N^ 1 is compared to the lower value of 
iV° 9M and avalanches occur until the new value for the 
critical parameter is reached and SOC is installed and 
maintained. 



SUMMARY AND OPEN ISSUES 



In all cases SOC is reached in a time period that de- 
pends on how close Ni is to the critical value. The num- 
ber of events needed to relax the system changes. 

A more realistic simulation has to account for when the 
superradiance condition is no longer fulfilled. It is then 
we expect that SOC to cease, followed by a transitory 
period, and afterwards the system will again settle into 
a SOC state characterized by the new parameters. Due 
to the limited computational capabilities and the very 
slow growth rate, we use the fact that the system will 
eventually be in SOC and initialize the automaton very 
close to its critical parameter, x — 0.99. 

To compute how much mass the BH has lost because of 
superradiance, we assume as an example that superradi- 
ance stops when 10% of the BH mass is lost, i.e. superra- 
diance stops for Mbh = M but immediately starts again 



Starting from rigorous analytical results describing the 
dynamics of a boson-cloud-black-hole (BC-BH) system, 
we derived a simple map from this dynamics to a cellular 
automaton. Within our working framework (see assump- 
tions at the end of Section [Hi] and Table HJ we have shown 
that the BC-BH system reaches self-organized criticality 
(SOC) and stays there as long as the superradiance con- 
dition is fulfilled. Since laboratory observations of the 
growth and collapse of Bose-Einstein condensates lead to 
a similar picture — albeit only for short timescales, see 
for instance Fig. 3 in [24j — we have some confidence that 
this is a physical property of our model and not an ar- 
tifact of inadequate modelling. When the superradiance 
condition is no longer valid, the evolution rules change 
quantitatively and the simulated system exhibits a tran- 
sition after which it again settles into SOC. By varying 
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FIG. 4: The distribution of event size, D(s) for a QCD-axion 
cloud in superradiance around a BH. The cloud was initialized 
at 99% of its critical boson number (the constraint Eq. (|17|l is 
obeyed). Transition to another set of superradiance param- 
eters is visible. The flow at criticality is given by rules (|19[) 
(upper panel) and for rules (|20p (lower panel). 



the rules we checked that this is a rule- independent result 
(in the sense that it is recovered for all the sets of rules 
that we have employed) which gives some confidence that 
this result holds for actual axion cloud near black holes. 
The quantitative details of the transition are rule de- 
pendent (Figs. [2J3]). For the first set of transition rules, 
Eq. ([l~9|) . the temporal width of the transition period is 
quite large, 25 x 10 4 r 9 . For the second set of transition 
rules, Eq. (f2"U)) . the temporal width of the transition is of 



the order 10 



The ratio of the medium numerical value 



of D(s) before and after transition is approximately 2 for 
four different evolutions of the system at criticality. The 
open issue is whether or not these aspects are related to 
underlying physics or are numerical artifacts. 

We assumed that a new superradiance domain is in- 
stantly valid, which led to very large number of events as- 
sociated with the transition when compared to the SOC 
state (Fig. 2]). An open question is whether or not this 



assumption is physically justified. This should be inves- 
tigated in future work. 

We expect that observational data of gravitational 
waves from a large range of masses for BH will be able 
to shed light on at least two very important issues. One 
is the physics at transition and observational data will 
help devise more realistic rules for suitable cellular au- 
tomata. Also, an appropriate observational database will 
help in determining the scaling function associated with 
this phenomena. 
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Appendix A: Sandpile model 

SOC may be studied through a simple sand-pile model 
(Fig. [5J in a CA framework [19(. This arrangement may 
be thought of as half of a symmetric sand pile with both 
ends open. The numbers z n characterizing the automa- 
ton represent height differences between successive posi- 
tions along the sand pile, z n = h(n) — h(n + 1). When 
placing a grain of sand in cell n the dynamics of this 
system follows the rules 



Z n r Z n -\- 1 

Z„_l — » Z n _l — 1. 



(Al) 



This model is a cellular automaton where the state of 
the discrete variable z n at time t + 1 depends on the state 
of the variable and its neighbors at time t. 

When the height difference becomes higher than a 
threshold (critical) value, z c , one unit of sand tumbles 
to the lower level 



Zn 
Zn±l 



z n z 

Zn±l + 1 for 2t, 



(A2) 



> z c 



The process continues until all the z n are below z c . At 
this point another grain of sand is added at a random 
site n through Eq. (jAll) . 
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FIG. 5: One dimensional sand pile automaton [19 |. 



Starting from a set of initial conditions, the pile grows 
and in the mean time the slope of the pile increases. It is 
the value of this slope, as seen in the Eqs. (|Aip . (IA2p . that 
will at some point reach a critical value. The meaning of 
the critical value is that if this value is reached and one 
more grain is added, material will slide. In the critical 



state the system is just barely stable with respect to fur- 
ther perturbations. Information about the perturbation 
will be known across all lengthscales in the system. 

An event consists of the seed (i.e. placing of the 
grain) and the subsequent information (grain) transfer 
or avalanche until complete relaxation. Each event has 
an associated time T (i.e. the time it takes the perturba- 
tion to die out), a cluster size s (i.e. the total number of 
slidings induced by the perturbation) and a total energy 
(i.e. release of energy, as transfers occur based on first 
principles that the final state is an energetically lower 
state). 

For this simple model, the distribution of event sizes 
D(s) is given by Q 



1 



D(s) = -e(L-s), 



(A3) 



where L is the system size and 0(x) is the Heaviside step- 
function. To put this in the form of Eq. ((TJ) , we rewrite 
it as 



D W = - 



H-i 



(A4) 



and by identification a = 1, t = 1, b = 1, a = 1 and the 

scaling function is Q{x) = x6 (1 — x). 
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